
clc;
clear all;
close all;

load ProgramImpact;

[row, col] = size(ProgramImpact);
Newmatrix(row, 24) = 0;


rou = 0.603; % estimated coefficient of lagged releases from STATA

col_year = 1;
col_part = 41;
col_facilityID = 48;
col_intrafirmspillover = 72;

col_prg = 9;   % level of releases
col_partmodified = 51;

col_totalimpact = 117;




index_year1990 = find(ProgramImpact(: ,col_year) == 1990);
tempMatrix1990 = ProgramImpact(index_year1990, :);
facilityIDs = unique(tempMatrix1990(:,col_facilityID));   


base = 1;
for i = 1: length(facilityIDs)
    facilityID = facilityIDs(i);
    Newmatrix(base, 1) = facilityID;   % column 1: facilityID
    index_facilityID = find(ProgramImpact(: ,col_facilityID) == facilityID);  % find across the whole sample
    tempMatrix_facilityID = ProgramImpact(index_facilityID, :);   % the matrix of this specific facilityID for all years
    
    indexpartYears = find((tempMatrix_facilityID(:, col_part) == 1) & (tempMatrix_facilityID(:, col_year) <=1995) );
    if (~isempty(indexpartYears))
       joinyear = min(tempMatrix_facilityID(indexpartYears, col_year));
       Newmatrix(base, 2) = joinyear;
    else
       Newmatrix(base, 2) = 0;  % 0 means never joined
    end
    
    indexintraspilloverYears = find((tempMatrix_facilityID(:, col_intrafirmspillover)>0) & (tempMatrix_facilityID(:, col_year)<=1996) );   % 1992 means spillover from 1991
    if (~isempty(indexintraspilloverYears))
       indexintraspilloverStartyear = min(tempMatrix_facilityID(indexintraspilloverYears, col_year));
       Newmatrix(base, 3) = indexintraspilloverStartyear;   % 1992 means spillover from 1991
    else
       Newmatrix(base, 3) = 0;  
    end   
    
    % 1990 releases
    release1990 = tempMatrix_facilityID( find(tempMatrix_facilityID(:,col_year) == 1990), col_prg);  % actual releases 1990
    Newmatrix(base, 4) = release1990;  
    
    index1991 =  find(tempMatrix_facilityID(:,col_year) == 1991);
    if (~isempty(index1991))
         impact1991 = tempMatrix_facilityID(index1991, col_totalimpact);  
         delta1991 = release1990* impact1991;
     else
         impact1991 = 0;
         delta1991= 0;
    end
    release1991 = release1990 + delta1991;
   
    
    index1992 =  find(tempMatrix_facilityID(:,col_year) == 1992);
    if (~isempty(index1992))
         impact1992 = tempMatrix_facilityID(index1992, col_totalimpact);   
         delta1992 = release1991* (impact1992+ impact1991*rou^(1992-1991));
     else
         impact1992 = 0;
         delta1992 = 0;
    end
    release1992 = release1990 + delta1991 + delta1992;
    
    
    
    index1993 =  find(tempMatrix_facilityID(:,col_year) == 1993);
    if (~isempty(index1993))
         impact1993 = tempMatrix_facilityID(index1993, col_totalimpact);   
         delta1993 = release1992* (impact1993+ impact1992*rou^(1993-1991));
     else
         impact1993 = 0;
         delta1993 = 0;
    end
    release1993 = release1990 + delta1991+ delta1992 + delta1993;
   
    
    
    index1994 =  find(tempMatrix_facilityID(:,col_year) == 1994);
    if (~isempty(index1994))
         impact1994 = tempMatrix_facilityID(index1994, col_totalimpact);   
         delta1994 = release1993* (impact1994+ impact1993*rou^(1994-1991));
     else
         impact1994 = 0;
         delta1994 = 0;
    end
    release1994 = release1990 + delta1991 + delta1992 + delta1993 + delta1994;
    
    
    
    index1995 =  find(tempMatrix_facilityID(:,col_year) == 1995);
    if (~isempty(index1995))
         impact1995 = tempMatrix_facilityID(index1995, col_totalimpact);   %v1992
         delta1995 = release1994* (impact1995+ impact1994*rou^(1995-1991));
     else
         impact1995 = 0;
         delta1995 = 0;
    end
    release1995 = release1990 + delta1991+ delta1992+ delta1993+ delta1994+ delta1995;
    
    
    index1996 =  find(tempMatrix_facilityID(:,col_year) == 1996);
     if (~isempty(index1996))
         impact1996 = tempMatrix_facilityID(index1996, col_totalimpact);   %v1992
         delta1996 = release1995* (impact1996+ impact1995*rou^(1996-1991));
     else
         impact1996 = 0;
         delta1996 = 0;
     end
    release1996 = release1990 + delta1991+ delta1992+ delta1993+ delta1994+ delta1995 + delta1996;
    
    
    Newmatrix(base, 5) = delta1991;
    Newmatrix(base, 6) = release1991;
    Newmatrix(base, 7) = delta1992;
    Newmatrix(base, 8) = release1992;
    Newmatrix(base, 9) = delta1993;
    Newmatrix(base, 10) = release1993;
    Newmatrix(base, 11) = delta1994;
    Newmatrix(base, 12) = release1994;
    Newmatrix(base, 13) = delta1995;
    Newmatrix(base, 14) = release1995;
    Newmatrix(base, 15) = delta1996;
    Newmatrix(base, 16) = release1996;
    if (Newmatrix(base, 4) >0)    % release1990;  
        Newmatrix(base, 17) = 100* delta1991 / Newmatrix(base, 4);
        Newmatrix(base, 18) = 100* delta1992 / Newmatrix(base, 4);
        Newmatrix(base, 19) = 100* delta1993 / Newmatrix(base, 4);
        Newmatrix(base, 20) = 100* delta1994 / Newmatrix(base, 4);
        Newmatrix(base, 21) = 100* delta1995 / Newmatrix(base, 4);
        Newmatrix(base, 22) = 100* delta1996 / Newmatrix(base, 4);
        Newmatrix(base, 23) = 100* (delta1991+ delta1992+ delta1993+ delta1994+ delta1995) / Newmatrix(base, 4);
        Newmatrix(base, 24) = 100* (delta1991+ delta1992+ delta1993+ delta1994+ delta1995+ delta1996) / Newmatrix(base, 4);
    end        
    
   base = base +1;
   
end  % for i = 1: length(facilityIDs)
    
Newmatrix(base: row, :) = [];

ImpactEachFacility = Newmatrix;
save ImpactEachFacility.mat ImpactEachFacility;
xlswrite('ImpactEachFacility',ImpactEachFacility);    
 


clear index_everjoin;

% for participants 
Impactmatrix_part = zeros(19, 7);
Impactmatrix_part(4: 15,1) = [1991 1991 1992 1992 1993 1993 1994 1994 1995 1995 1996 1996]';
%%%% for those ever joined  %%%%%%%
index_everjoin = find(Newmatrix(:,2)>0);  % col 2 is joined year
obsN = length(index_everjoin);   % observation number 
Impactmatrix_part(1,2) = obsN;
indexNonzon = (Newmatrix(index_everjoin, 4)>0); % only sum over those releases>0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 4))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(2,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 4));   % 1990 average releases
Impactmatrix_part(3,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 4))*(length(indexCombine))^(-0.5);  % s.e.
end

% reduction in each year
indexNonzon = (Newmatrix(index_everjoin, 17)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 17))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(4,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 17)); 
Impactmatrix_part(5,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 17))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 18)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 18))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(6,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 18)); 
Impactmatrix_part(7,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 18))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 19)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 19))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(8,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 19)); 
Impactmatrix_part(9,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 19))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 20)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 20))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(10,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 20)); 
Impactmatrix_part(11,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 20))*(length(indexCombine))^(-0.5);  % s.e.
end


indexNonzon = (Newmatrix(index_everjoin, 21)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 21))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(12,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 21)); 
Impactmatrix_part(13,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 21))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 22)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 22))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(14,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 22)); 
Impactmatrix_part(15,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 22))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 23)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 23))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(16,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 23)); 
Impactmatrix_part(17,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 23))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 24)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 24))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_part(18,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 24)); 
Impactmatrix_part(19,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 24))*(length(indexCombine))^(-0.5);  % s.e.
end



%%%% do separately for joined in different years
for j = 1991: 1995
    index_everjoin = find(Newmatrix(:,2) == j);  
    obsN = length(index_everjoin);   % observation number 
    Impactmatrix_part(1,2+ (j-1990)) = obsN;
    indexNonzon = (Newmatrix(index_everjoin, 4)>0); % only sum over those releases>0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 4))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(2,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 4));   % 1990 average releases
    Impactmatrix_part(3,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 4))*(length(indexCombine))^(-0.5);  % s.e.
    end
    
    % reduction in each year
    indexNonzon = (Newmatrix(index_everjoin, 17)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 17))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(4,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 17)); 
    Impactmatrix_part(5,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 17))*(length(indexCombine))^(-0.5);  % s.e.
    end

    indexNonzon = (Newmatrix(index_everjoin, 18)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 18))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(6,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 18)); 
    Impactmatrix_part(7,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 18))*(length(indexCombine))^(-0.5);  % s.e.
    end
    
    indexNonzon = (Newmatrix(index_everjoin, 19)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 19))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(8,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 19)); 
    Impactmatrix_part(9,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 19))*(length(indexCombine))^(-0.5);  % s.e.
    end


    indexNonzon = (Newmatrix(index_everjoin, 20)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 20))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(10,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 20)); 
    Impactmatrix_part(11,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 20))*(length(indexCombine))^(-0.5);  % s.e.
    end


    indexNonzon = (Newmatrix(index_everjoin, 21)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 21))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(12,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 21)); 
    Impactmatrix_part(13,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 21))*(length(indexCombine))^(-0.5);  % s.e.
    end
    
    indexNonzon = (Newmatrix(index_everjoin, 22)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 22))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(14,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 22)); 
    Impactmatrix_part(15,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 22))*(length(indexCombine))^(-0.5);  % s.e.
    end

    indexNonzon = (Newmatrix(index_everjoin, 23)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 23))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(16,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 23)); 
    Impactmatrix_part(17,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 23))*(length(indexCombine))^(-0.5);  % s.e.
    end


    indexNonzon = (Newmatrix(index_everjoin, 24)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 24))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_part(18,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 24)); 
    Impactmatrix_part(19,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 24))*(length(indexCombine))^(-0.5);  % s.e.
    end
    
end




%%% for non-participants face intra-firm spillover
Impactmatrix_intraspillover = zeros(19, 7);
Impactmatrix_intraspillover(4: 15,1) = [1991 1991 1992 1992 1993 1993 1994 1994 1995 1995 1996 1996]';
%%% for those ever joined  %%%%%%%
index_everjoin = find(Newmatrix(:,3)>0);   % col 3 is intra-spillover start year!!
obsN = length(index_everjoin);   % observation number 
Impactmatrix_intraspillover(1,2) = obsN;
indexNonzon = (Newmatrix(index_everjoin, 4)>0); % only sum over those releases>0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 4))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(2,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 4));   % 1990 average releases
Impactmatrix_intraspillover(3,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 4))*(length(indexCombine))^(-0.5);  % s.e.
end
    
% reduction in each year
indexNonzon = (Newmatrix(index_everjoin, 17)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 17))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(4,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 17)); 
Impactmatrix_intraspillover(5,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 17))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 18)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 18))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(6,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 18)); 
Impactmatrix_intraspillover(7,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 18))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 19)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 19))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(8,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 19)); 
Impactmatrix_intraspillover(9,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 19))*(length(indexCombine))^(-0.5);  % s.e.
end


indexNonzon = (Newmatrix(index_everjoin, 20)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 20))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(10,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 20)); 
Impactmatrix_intraspillover(11,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 20))*(length(indexCombine))^(-0.5);  % s.e.
end



indexNonzon = (Newmatrix(index_everjoin, 21)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 21))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(12,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 21)); 
Impactmatrix_intraspillover(13,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 21))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 22)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 22))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(14,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 22)); 
Impactmatrix_intraspillover(15,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 22))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_everjoin, 23)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 23))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(16,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 23)); 
Impactmatrix_intraspillover(17,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 23))*(length(indexCombine))^(-0.5);  % s.e.
end


indexNonzon = (Newmatrix(index_everjoin, 24)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_everjoin, 24))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_intraspillover(18,2) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 24)); 
Impactmatrix_intraspillover(19,2) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 24))*(length(indexCombine))^(-0.5);  % s.e.
end




%%% do separately for joined in different years
for j = 1991: 1995
    index_everjoin = find(Newmatrix(:,3) == j+1);     % column 3 is intra-spillover year
    obsN = length(index_everjoin);   % observation number 
    Impactmatrix_intraspillover(1,2+ (j-1990)) = obsN;
    indexNonzon = (Newmatrix(index_everjoin, 4)>0); % only sum over those releases>0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 4))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(2,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 4));   % 1990 average releases
    Impactmatrix_intraspillover(3,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 4))*(length(indexCombine))^(-0.5);  % s.e.
    end

%    reduction in each year
    indexNonzon = (Newmatrix(index_everjoin, 17)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 17))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(4,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 17)); 
    Impactmatrix_intraspillover(5,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 17))*(length(indexCombine))^(-0.5);  % s.e.
    end

    indexNonzon = (Newmatrix(index_everjoin, 18)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 18))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(6,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 18)); 
    Impactmatrix_intraspillover(7,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 18))*(length(indexCombine))^(-0.5);  % s.e.
    end

    indexNonzon = (Newmatrix(index_everjoin, 19)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 19))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(8,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 19)); 
    Impactmatrix_intraspillover(9,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 19))*(length(indexCombine))^(-0.5);  % s.e.
    end


    indexNonzon = (Newmatrix(index_everjoin, 20)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 20))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(10,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 20)); 
    Impactmatrix_intraspillover(11,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 20))*(length(indexCombine))^(-0.5);  % s.e.
    end


    indexNonzon = (Newmatrix(index_everjoin, 21)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 21))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(12,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 21)); 
    Impactmatrix_intraspillover(13,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 21))*(length(indexCombine))^(-0.5);  % s.e.
    end

    indexNonzon = (Newmatrix(index_everjoin, 22)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 22))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(14,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 22)); 
    Impactmatrix_intraspillover(15,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 22))*(length(indexCombine))^(-0.5);  % s.e.
    end
    
    indexNonzon = (Newmatrix(index_everjoin, 23)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 23))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(16,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 23)); 
    Impactmatrix_intraspillover(17,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 23))*(length(indexCombine))^(-0.5);  % s.e.
    end

    indexNonzon = (Newmatrix(index_everjoin, 24)<0); % only sum over those reduction < 0
    indexNonNAN = (~isnan(Newmatrix(index_everjoin, 24))==1);
    indexCombine = (indexNonzon & indexNonNAN);
    if (~isempty(indexCombine))
    Impactmatrix_intraspillover(18,2+ (j-1990)) = mean(Newmatrix(index_everjoin(find(indexCombine==1)), 24)); 
    Impactmatrix_intraspillover(19,2+ (j-1990)) = std(Newmatrix(index_everjoin(find(indexCombine==1)), 24))*(length(indexCombine))^(-0.5);  % s.e.
    end
    
end



% for non-participants face intra- and inter-firm spillover
Impactmatrix_interspillover = zeros(19, 2);
Impactmatrix_interspillover(4: 15,1) = [1991 1991 1992 1992 1993 1993 1994 1994 1995 1995 1996 1996]';
%%%% for those ever joined  %%%%%%%
%index_interspillover = find( (Newmatrix(:,2)==0) & (Newmatrix(:,3)==0) ); 
index_interspillover = find( (Newmatrix(:,2)==0) & (Newmatrix(:,3)>=0) ); 
obsN = length(index_interspillover);   % observation number 
Impactmatrix_interspillover(1,2) = obsN;
indexNonzon = (Newmatrix(index_interspillover, 4)>0); % only sum over those releases>0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 4))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(2,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 4));   % 1990 average releases
Impactmatrix_interspillover(3,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 4))*(length(indexCombine))^(-0.5);  % s.e.
end

% reduction in each year
indexNonzon = (Newmatrix(index_interspillover, 17)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 17))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(4,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 17)); 
Impactmatrix_interspillover(5,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 17))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_interspillover, 18)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 18))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(6,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 18)); 
Impactmatrix_interspillover(7,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 18))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_interspillover, 19)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 19))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(8,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 19)); 
Impactmatrix_interspillover(9,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 19))*(length(indexCombine))^(-0.5);  % s.e.
end


indexNonzon = (Newmatrix(index_interspillover, 20)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 20))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(10,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 20)); 
Impactmatrix_interspillover(11,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 20))*(length(indexCombine))^(-0.5);  % s.e.
end



indexNonzon = (Newmatrix(index_interspillover, 21)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 21))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(12,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 21)); 
Impactmatrix_interspillover(13,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 21))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_interspillover, 22)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 22))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(14,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 22)); 
Impactmatrix_interspillover(15,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 22))*(length(indexCombine))^(-0.5);  % s.e.
end

indexNonzon = (Newmatrix(index_interspillover, 23)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 23))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(16,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 23)); 
Impactmatrix_interspillover(17,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 23))*(length(indexCombine))^(-0.5);  % s.e.
end


indexNonzon = (Newmatrix(index_interspillover, 24)<0); % only sum over those reduction < 0
indexNonNAN = (~isnan(Newmatrix(index_interspillover, 24))==1);
indexCombine = (indexNonzon & indexNonNAN);
if (~isempty(indexCombine))
Impactmatrix_interspillover(18,2) = mean(Newmatrix(index_interspillover(find(indexCombine==1)), 24)); 
Impactmatrix_interspillover(19,2) = std(Newmatrix(index_interspillover(find(indexCombine==1)), 24))*(length(indexCombine))^(-0.5);  % s.e.
end




ReductionPercentage = zeros(57, 7);
ReductionPercentage(1: 19, 1: 7) = Impactmatrix_part;
%ReductionPercentage(20: 38, 1: 7) = Impactmatrix_intraspillover;  % ignore these rows
ReductionPercentage(39: 57, 1: 2) = Impactmatrix_interspillover;  % refer to inter and intra-firm spillover



save ReductionPercentage.mat ReductionPercentage;
xlswrite('ReductionPercentage',ReductionPercentage);    
        

